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ABSTRACT 

The  Future  Combat  Systems  (FCS)  will  transform 
the  Army  into  a  faster,  more  agile  force  with  more 
lethal  and  more  survivable  capabilities.  The  FCS  re¬ 
quires  computational  technologies  that  can  simulate 
and  assess  lethality  and  survivability  of  the  combat 
systems  under  real  battlefield  situations  involving  the 
global  structural  dynamics  at  system  level,  and  highly 
localized  contact-impact,  damage,  failure,  penetration, 
and  fragmentation.  In  order  to  integrate  the  aforemen¬ 
tioned  capabilities  into  a  standalone  computational 
code,  a  computational  framework  taking  into  consid¬ 
eration  different  mechanics  descriptions  and  formula¬ 
tions,  different  physical  events  such  as  structural  dy¬ 
namics,  contact- impact,  damage,  failure,  penetration, 
and  fragmentation  need  to  be  carefully  addressed.  In 
this  regard,  different  computational  formulations  in 
space  such  as  finite  element  methods  and  meshless 
methods,  and  in  time  such  as  implicit  algorithms,  ex¬ 
plicit  algorithms,  and  the  like  need  to  be  seamlessly 
integrated.  Existing  technologies  focus  on  an  isolated 
physical  event  or  a  limited  coupled  physical  event.  A 
computational  framework  that  can  integrate  the  afore¬ 
mentioned  multi-physical  events  for  combat  systems 
commonly  encountered  in  real  battlefield  situations 
does  not  exist  to-date.  We  address  the  above  issues 
with  applications  to  engineering  problems,  and  the 
present  approach  compares  well  with  experimental  re¬ 
sults. 

1.  INTRODUCTION 

Hydrocodes  are  important  tools  for  survivability  and 
lethality  of  armor/anti-armor  applications  [Schraml 
et  ah,  2002].  The  existing  hydrocodes  mostly  utilize 
fracture  mechanics  based  failure  models  (the  failure 
of  the  material  is  determined  by  the  maximum  accu¬ 
mulated  plastic  strain)  [Zukas,  2004].  In  contrast  to 


the  fracture  mechanics  based  failure  models,  which  fo¬ 
cus  on  when  the  material  fails,  the  continuum  dam¬ 
age  mechanics  based  failure  models  (the  failure  of  the 
material  is  determined  by  the  evolution  of  an  inter¬ 
nal  variable  that  measures  the  micro-defect  density) 
can  further  provide  an  in-depth  understanding  of  the 
physical  mechanisms  of  how  the  material  fails  [Krajci- 
novic,  2000].  Furthermore,  since  Lemaitre  introduced 
an  approach  to  incorporate  fracture  mechanics  based 
failure  model  into  the  continuum  damage  mechan¬ 
ics  based  isotropic  damage  failure  models  [Lemaitre, 
1985a,  Lemaitre,  1985b],  the  Lemaitre  damage  model 
has  been  successfully  applied  for  the  quasi-static  fail¬ 
ure  of  steel  [Celentano  et  ah,  2004]  and  high-velocity 
penetration  [Bprvik  et  ah,  2001].  The  present  study 
aims  to  extend  the  Lemaitre  damage  model  to  the  hy¬ 
drodynamic  type  of  constitutive  model  for  high  veloc¬ 
ity  penetration  problems. 

With  the  objective  towards  the  development  of  a 
computational  framework  that  can  simulate  the  inte¬ 
grated  multi-physics  and  multi-time  scale  battlefield 
events,  the  development  of  highly  efficient  and  scalable 
multi-time  scale  and  adaptive  semi-implicit /explicit 
hybrid  dynamic  computational  is  the  key  for  predict¬ 
ing  both  the  long  term  global  dynamic  response  of 
structural  systems  such  as  the  entire  transient  response 
of  the  combat  system  and  structural  configuration, 
in  conjunction  with  effects  imparted  due  to  localized 
short  term  dynamic  damage  due  to  threats  such  as 
penetration.  Quite  often,  a  representative  section  of 
the  vehicle  component  is  mostly  considered  (for  ex¬ 
ample  a  small  armor  plate  section)  for  the  simulation 
in  survivability /lethality  designs.  Increasing  model  fi¬ 
delity  requires  high  mesh  resolutions.  Longer  simula¬ 
tion  times  are  additionally  needed  to  assess  post  dy¬ 
namic  simulations  at  the  system  level  due  to  localized 
contact-impact-penetration-damage  events.  Compu¬ 
tational  requirements  grow  dramatically  as  the  size 
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and  complexity  of  the  computational  models  increase. 
As  the  element  size  is  decreased  to  capture  the  finer 
details  such  as  large  deformations  in  geometry  and  for 
overall  response,  the  time  step  also  significantly  de¬ 
creases.  This  is  especially  cumbersome  and  also  im¬ 
practical  for  large-scale  simulations  such  as  an  entire 
Army  vehicle  system  to  be  completed  in  reasonable 
time.  Using  traditional  explicit  and  implicit  finite  el¬ 
ement  codes,  the  coupling  of  standalone  explicit  and 
implicit  codes  is  not  feasible  and/or  practical  due  to: 
1)  approximation  errors  in  transferring  results  from 
one  code  to  the  other,  and  2)  impractical  to  do  data 
transfer  if  damage,  penetration  and  fragmentation  are 
involved  with  flying  debris.  Towards  this  end,  the 
present  efforts  circumvent  several  existing  drawbacks 
via  new  computational  methods  and  algorithms  and 
include  the  total  integration  and  fusion  of  various  com¬ 
putational  technologies.  Therefore,  an  approach  of 
adaptive  multi-scale  semi-implicit /explicit  hybrid  for¬ 
mulations  and  computational  algorithms  would  capi¬ 
talize  upon  both  the  advantages  of  implicit  and  explicit 
time  integration  algorithms  to  avoid  the  existing  draw¬ 
backs  and  is  desirable  for  the  simulations  of  numerous 
applications  involving  lethality  and  survivability.  Such 
a  hybrid  approach  is  not  feasible  within  the  standard 
framework  of  linear  multi-step  (LMS)  methods  that 
exist  in  most  Lagrangian  codes  [Zhou  and  Tamma, 
2004].  The  approach  for  developing  hybrid  and  adap¬ 
tive  multi-scale  computational  algorithms  is  based  on 
the  recent  and  highly  successful  developments  of  ro¬ 
bust  element  based  smart-switch  algorithms  in  time 
which  smartly  switch  on/off  the  necessary  algorithms 
on  an  element  basis  for  the  given  simulation  model 
(between  a  semi-implicit  formulation  which  do  not  re¬ 
quire  any  nonlinear  iterations  [Zhou  et  ah,  2005]  and 
a  explicit  time  formulations  [Sha  et  ah,  1996],  both 
of  which  have  been  derived  from  a  novel  theory  on 
the  time  dimension).  The  automated  built-in  switch 
is  based  on  the  element  physical  and  geometrical  mea¬ 
sures,  and  on  the  state  of  the  damage  variables;  as 
a  consequence,  if  the  time  steps  for  explicit  elements 
in  the  numerical  mesh  become  prohibitively  small  due 
to  severe  large  deformations  and  distortions,  only  the 
corresponding  elements  are  switched  to  semi-implicit 
formulations  with  a  much  larger  time  step.  Time  adap¬ 
tive  features  are  embedded  into  the  framework  to  pro¬ 
vide  optimal  time  steps  for  assessing  both  short/long 
term  dynamic  simulations.  Since  the  total  simulation 
framework  of  semi-implicit  time  formulations  is  read¬ 
ily  integrated  to  the  naturally  scalable  explicit  formu¬ 
lations  while  the  load  balancing  is  performed  appropri¬ 
ately,  the  overall  HPC  simulations  can  readily  handle 
the  multiple  time  scale  applications. 

The  presentation  of  this  exposition  is  outlined  as 
the  follows.  First  the  dynamic  equilibrium  equation  in¬ 


cluding  frictional  contact  boundary  conditions  in  Eule- 
rian  form  is  formulated.  Then  the  deviatoric  portion  of 
the  hydrodynamic  constitutive  equation  incorporating 
the  elasto-visco-thermo-plastic-damage  is  formulated. 
And,  the  hydrostatic  portion  of  the  hydrodynamic  con¬ 
stitutive  equation  utilizing  the  equation  of  state  and 
internal  energy  equation  is  then  discussed.  The  time 
integration  procedure  is  described.  Finally,  numeri¬ 
cal  validation  are  performed  to  demonstrate  the  appli¬ 
cability  and  robustness  of  the  present  computational 
framework  and  formulation. 

2.  EQUATION  OF  MOTION 

Let  the  open  set  Qt  C  5ft3  be  the  domain  of  inter¬ 
est  of  the  configuration  at  time  £,  with  the  bound¬ 
ary  dflt  =  r{  U  r/,  r{  n  Tf  =  0,  and  the  closure 
Qt  =  Qt  U  dQf  For  finite  strain  dynamic  problems, 
the  dynamic  equilibrium  equation  of  the  configuration 
at  time  t  in  the  Eulerian  form  with  Cauchy  stress  is 
described  by 

bt  in  Qt  (1) 

/  on  r{  (2) 
xd  on  Tf  (3) 

where  pt  is  the  material  density  of  the  configuration 
at  time  £,  vt  is  the  material  particle  velocity  at  time  £, 
at  is  the  Cauchy  stress  of  the  configuration  at  time  £, 
rj  is  the  viscous  damping  ratio  in  the  sense  of  Rayleigh 
damping,  bt  is  the  body  force,  n  is  the  boundary  sur¬ 
face  outward  unit  vector,  xt  is  material  particle  coordi¬ 
nate,  r{  is  the  boundary  subject  to  traction  boundary 
condition,  Yd  is  the  boundary  subject  to  prescribed 
displacement  boundary  condition. 

Let  an  =  an  =  {crnk}  =  -{^(7^^}  G  5ft3  be 
the  normal  stress  on  the  contact  surface,  n  be  the  nor¬ 
mal  direction  vector  on  the  contact  surface,  g  G  5ft3  be 
the  relative  gap  vector  between  point  pairs  on  the  con¬ 
tact  surfaces.  Let  gn  and  gT  be  the  normal  and  tan¬ 
gential  directions  of  the  gap  vector  g ,  then  g  =  gn+gT 
and  gn  gT  =  0.  And,  the  impenetrable  constraint  in 
the  normal  direction  on  the  contact  surfaces  is  given 
as  the  following  Kuhn- Tucker  complementarity  condi- 
tions,  (t n, gn  -  \cnn)Tc{t)  =  0,  \cn  >  0,  = 

<T  ■  n  >  0,  A cni>n{(Tn)  =  0,  V  Tn  e  1pn{rn)  >  o, 
where  >  0  represents  the  impenetrability  condition 
on  the  contact  surface,  and  (•,  »)r  =  fr  •  :  •ds  repre¬ 
sents  the  inner-product. 

Let  n  =  {V  |  V  G  5ft3  x  [0,  oo)},  then  a  yield 
function  of  the  frictional  stress  vector  in  the  tangen¬ 
tial  direction  on  the  contact  surface  can  be  defined  by 
2pr(ar)  =  gc  |  an  \  —  \  aT  |>  0  on  Tc  x  [0,  oo),  where 
<xT  =  {crtfc}  =  -{niaik  -  niakjUjHk}  €  SR3  is  the  tan- 
gential  stress  vector  on  the  contact  surface,  and  fic  G  5ft 
is  the  Coulomb  frictional  coefficient. 


d(ptvt) 

dt 


■  ptTJVt  -  V  •  <Xt  = 

crtn  = 

xt  = 


2 


Let  gT  G  5ft3  be  the  the  relative  tangential  ve¬ 
locity  of  the  point  pairs  on  the  contact  surfaces, 
the  classical  Coulomb  friction  law  thus  is  repre¬ 
sented  by  (rr  —  crT,gT)  <  0,  crr  G  2pr(ar)  >  0, 
V  rr  G  Vv(rr)  >  0.  The  corresponding  weak 
form  employing  the  Lagrangian  multiplier  is  given 
as  (rT,gT  -  A cTVaTipT(crT))rc(t)  =  0,  ipT(crT)  >  0, 
\(T  >  0,  \^iJjT(crT)  =  0,  V  rT  G  ipT(rT)  >  0,  where 

Define  the  following  set  of  virtual  displacement 
field  as,  [iLo(Dt)]3  :=  {w  \  w  G  [iL1(^)]3,  w  = 
0  on  Td},  and  the  set  of  virtual  velocity  field  as, 
:=  {v  |  v  G  [iJ1($7)]3,  v  =  xdt  on  Td}, 
where  H 1(Q)  is  the  first  order  differentiable  Hilbert 
space.  The  Euler ian  weak  form  of  the  dynamic  equi¬ 
librium  equation  subjected  to  the  contact  boundary 
condition  is  given  as  the  following. 

(w,  ■  +  pgv) nt  +  (T>(w),  a)Qt 

=  (w,b)nt+(w,f)rf{t)+(Sg,crn+aT)  rc(t)  (4) 
(Tn,9n-  Kn)rc(t)  =0  (5) 

(tt ,9t  -  A£WrVv(ov))rc(t)  =  o  (6) 

and  subject  to  the  constraints  A^  >  0,  2pn(crn)  = 
<Tn-n>  0,  \^!pn(crn)  =  0,  on  Tc(t),  and  i>T(crT)  >  0, 
K  >  0,  A cT^r((TT)  =  0,  on  Tc(t). 

In  the  matrix  form,  the  equation  of  motion  can  be 
written  as, 


Mii  +  Cu  +  Ku 

= 

f  (t)  +  gtst  +  gt 

N </>  (7) 

A 

= 

NtG(x  +  u)  >  0 

(8) 

> 

0 

(9) 

A  T  cf) 

0 

(10) 

st 

— -(g  T  +  AAt) 

(11) 

Ml <t>\a  -  IIstIU 

> 

0 

(12) 

where  M  is  the  mass  matrix,  C  is  the  damping  matrix, 
K  is  the  stiffness  matrix,  ii  is  the  acceleration,  u  is  the 
velocity,  and  u  is  the  displacement,  f  is  the  external 
loading,  G  is  the  constraint  matrix,  s t  is  the  tangen¬ 
tial  contact  force,  N  is  the  normal  gap  matrix,  <j>  is  the 
normal  contact  force  magnitude,  A  is  the  normal  gap 
magnitude,  x  is  the  position  vector,  fi  is  the  Coulomb 
coefficition  of  friction,  A  is  the  tangential  gap  matrix, 
and  e  is  an  infinitesimal  positive  number. 

3.  CONSTITUTIVE  EQUATION 


where  a  is  the  effective  stress  tensor,  a  is  the  Cauchy 
stress  tensor,  ip  is  the  isotropic  damage  state  variable, 
ajj  is  the  deviatoric  portion  of  the  effective  stress  ten¬ 
sor,  &h  is  the  hydrostatic  portion  of  the  effective  stress 
tensor,  G  is  the  shear  modulus,  is  the  elastic  part 
for  the  deviatoric  portion  of  the  velocity  strain  tensor, 
p  is  the  hydrostatic  pressure,  p  is  the  material  density, 
e  is  the  internal  energy  per  unit  mass,  I  is  the  rank  two 
identity  tensor,  and  the  symbol  «v  denotes  to  either 
the  Jaunamm,  the  Green-Naghdi,  or  Truesdell  stress 
rate. 

3.1  Damage  Models  for  the  Deviatoric  Stress 

From  a  practical  perspective,  the  Lemaitre  plastic- 
damage  model  [Lemaitre,  1996]  is  adapted  for  the 
deviatoric  stress  of  the  hydrodynamic  constitutive 
equation.  For  the  low-  and  high-velocity  im¬ 
pact/penetration  problems,  when  the  amount  of  dissi¬ 
pative  work  done  by  the  plasticity  and  damage  in  the 
system  does  not  raise  the  temperature  significantly, 
for  example,  a  projectile  impacting  on  a  thin  struc¬ 
ture,  the  Lemaitre  plastic-damage  model  should  be 
applicable.  However,  when  the  amount  of  dissipa¬ 
tive  work  gives  rise  to  a  rapid  increase  of  temperature 
for  the  material,  for  example,  in  the  case  of  hyperve¬ 
locity  penetration  problems,  a  temperature  dependent 
thermo- visco-plastic  model  such  as  the  Johnson-Cook 
fracture  model  [Johnson  and  Cook,  1985]  needs  to  be 
modified  and  incorporated  into  the  Lemaitre  plastic- 
damage  model  to  account  for  the  temperature  and 
strain  rate  effect.  The  modified  Johnson-Cook  thermo- 
visco-plastic-damage  model  is  described  next. 

The  von  Mises  tensile  flow  stress  of  the  Johnson- 
Cook  model  [Johnson  and  Cook,  1985]  is  given  by, 

ajC  =  (A  +  Ben)(l  +  Cln(-))(l-fm)  (14) 

where  A ,  H,  C,  n,  and  m  are  material  constants, 
£0  =  1.0s_1,  f  =  (T  -  Troorn) / (Tmeit  Troom),  Troorn 
is  the  room  temperature,  Tmeit  is  the  material  melting 
temperature.  And  the  fracture  strain  of  the  Johnson- 
Cook  model  [Johnson  and  Cook,  1985]  is  given  by 

eF  =  (D1+D2eD3^)(l  +  D4\n(-))(l  +  D5T)  (15) 
for  —  <  1.5,  where  D1-D5  are  material  constants. 

< 7  eg 

For  —  >1.5,  the  model  becomes 

&  eg 


Consider  the  isotropic  damage  hydrodynamic  consti¬ 
tutive  equation  given  by, 


u  =  - - =  CTd  +  CTh 

)  „  1 

|  y = 2 GV% 

{  aH  =  -p{p,e)I 


(13) 


eF  =  (D1  +  D2e1-5D*)(  1  +  D4  ln(>))(l  +  DbT)  (16) 

£o 


Let  the  thermo-visco-plastic-damage  potential 
function  be 


F'tvpdip' >  Rtypi^)  5  Q  (^)j  L,  (/?) 
=ftvP((T,  RtvP(e),  Q_1(e);  <p)  +  fd(Y,  <p) 
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where 


ftvp 

=  Oeq  -  Rtvp[e)  -  Q  1(i)  <  0 

(18) 

Rtvp 

=  (A  +  Ben)(l  -  Tm) 

(19) 

Q-+) 

II 

+ 

II 

O 

5 " 

S+v,. 

(20) 

frp 

Q(frp) 

=  e  =  i0eCRtvpH(frp) 

(21) 

Y  is  the  energy  release  rate  density,  fd  is  the  damage 
potential,  <p  is  the  damage  internal  variable,  and  H(m) 
is  the  Heaviside  function. 


Thus,  the  Modified  Johnson-Cook  thermo-visco- 
plastic-damage  model  yields  the  constitutive  equation 
in  consistent  form  as 


_v  _  f,  ,  („  C  :  ntvpd  0  ntvpd  :  C  \  ^ 

& D  v  I  ^  jjtvpd  .  £  .  rptvpd  _|_  pjtvpd  J  ’ 

(22) 


where 


HtvPd  =  h+c{h+^pp)Hi_) 
y  eA t  ’  +0 

H(e  -  eD)<?eq 


n 


(1  —  <p)2(£f  -  £d) 
h  =  Bnen-\l-fm) 

3cr  D 


tvp  _ 


2ae 


(23) 

(24) 

(25) 


Temperature  for  Adiabatic  Condition 

During  the  high  strain  rate  deformation,  the  tempera¬ 
ture  of  the  material  rises  due  to  the  dissipative  plastic 
work.  For  the  short  duration  of  the  material  deforma¬ 
tion  under  high  strain  rate,  the  thermal  equilibrium 
of  the  material  cannot  be  established.  Therefore,  the 
temperature  raise  can  be  computed  by  utilizing  the 
adiabatic  condition.  And,  the  temperature  equation  is 
given  as, 

pcT  =  pea  :  (VD  -  VeD)  (26) 

where  p  is  the  material  density,  c  is  the  material  heat 
capacity,  and  pe  =  0.9  is  efficiency  ratio  [Bammann 
et  ah,  1993]. 


rate  of  internal  energy  per  current  unit  mass  is  given 
as 

e  =  ~(crD  :  V -p  tr(T>))  (27) 

P 

where  p  is  the  hydrostatic  pressure,  and  tr(X>)  is  the 
trace  of  the  velocity  strain.  And,  the  Mie-Griineisen 
equation  of  state  [Walsh  et  ah,  1957]  considered  here 
is  given  by 


p(p ,  e)  =  (1  -  ^7 P)Ph  +  ipe  (28) 

where  p  =  p  —  1,  p  =  p/po,  Po  is  the  material  den¬ 
sity  at  the  initial  reference  state,  7  is  the  Gruneisen 
parameter,  70  is  the  Gruneisen  parameter  at  the  ini¬ 
tial  reference  state,  pn  is  the  Hugoniot  pressure  at  the 
density  p. 

Let  be  the  deviatoric  Cauchy  stress  at  time 

t  +  At  obtained  from  the  deviatoric  part  of  the  hydro- 
dynamic  constitutive  relation,  a  time  integration  al¬ 
gorithm  for  solving  the  hydrostatic  pressure  from  the 
equation  of  state  is  given  as  the  following. 


Algorithm  1 

For  the  time  increment  from  t  to  t+At,  the  hydrostatic 
pressure  can  be  integrated  from  following  algorithm. 


t+At  _ 


Leos  +  Beoset+At  pt+At 

1  +  3JBeos/?A7+At 


(29) 


where 

A 

s^-eos 

Beos 

~t+At 


(!  -  iph)VH 


=  7 


et  +  ■ 


1 


it+ At 


p 


(aertD  +  /3cr 


D 


(30) 

(31) 

*)  :  A£t+At 


-3ap*A£*+At] 

(32) 

/’t+At 

A£t+At  = 

J  Vdt ’ 

/■t+At 

(33) 

A£j+At  = 

J  Rndf 

(34) 

Vh  = 

-trace(T>) 

(35) 

crD  = 

(1  -  <p)crD 

(36) 

P  = 

1  —  a 

(37) 

3.2  Equation  of  State  for  the  Hydrostatic  Pres¬ 
sure 

The  equation  of  state  (EOS)  along  with  the  material 
internal  energy  equation  governs  the  relation  between 
the  pressure  and  density  implicitly,  and  serves  as  the 
constitutive  equations  which  express  the  state  that  the 
material  can  achieve  [Walsh  and  Christian,  1955].  The 


Artificial  Viscosity 

For  numerical  shock  wave  computation,  the  artificial 
viscosity  is  needed  to  stabilize  the  numerical  oscilla¬ 
tion  [von  Neumann  and  Richtmyer,  1950].  The  classi¬ 
cal  artificial  viscosity  is  adapted  as  the  following. 
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aq  Ax2pVj{  +  02  A xpVn  T>h  <  0 
0  VH  >  0 


(38) 


where  oq  and  <a2  are  the  artificial  viscosity  parameters, 
Ax  is  the  element  characteristic  length.  Therefore,  the 
internal  energy  rate  and  the  Cauchy  stress  is  modified 
as 

e  =  a  :  T>  -3(p  +  q)VH  (39) 

cr  =  crD-(p  +  q)I  (40) 

4.  TIME  INTEGRATION 

The  implementation  procedures  for  updated 
largrangian  formulation  of  nonlinear  dynamic  system 
with  frictional  contact  boundary  is  described  in  the 
following. 

Time  Integration  1 

Preparing  for  start:  set  n  =  0,  if  the  element  time 
step  is  greater  than  the  critical  time  step  of  central 
difference  method,  set  £  =  0.58383879542,  otherwise 
set  £  =  0,  and  let  matrix  Q  is  updated  for  every  j 
step;  Decompose  the  mass  matrix  M  by  M  =  LTL, 
and  let  q  =  Lu. 

1°.  If  n  =  0  or  n  is  integer  time  of  the  j  then  do  fol¬ 
lowing:  set  p  =  n  and  generate  the  tangential  stiffness 

* 

matrix  Kp,  form  the  matrix  Kp,  Qp,  C,  M  as 


Kp 

=  L-TKpL-1 

(41) 

Q  p 

- 

(42) 

c 

=  L  7  CL  1 

(43) 

* 

M 

=  I+^AiC 

(44) 

5°.  Find  0,  Aqn+2,  qn+i  and  qn+i  by  using  the 
linear  complementary  bi-projection  (LCBP)  method 
to  solve  complementary  equation: 

A(0)  =  NTG(xn+i  +  Aqn+2(0))  >  0  (52) 

0  >  0  (53) 

XT(<P)<P  =  0  (54) 

6°.  n  :=  n  +  1  goto  2°. 

Remark  0.1 

1.  The  algorithm  is  L-stable  and  second  or¬ 
der  accuracy  for  nonlinear  dynamic  system 
with  frictional  contact  boundary  when  £*  > 
0.583838795422. 

2.  The  algorithm  is  the  central  difference  method 
when  £  =  0. 

5.  Numerical  Validation 

The  numerical  example  considers  experiments  con¬ 
ducted  by  Forrestal  and  Hanchak  with  blunt  4340  Rc 
38  steel  cylinder  projectiles  impacting  circular  HY-100 
steel  plates  [Forrestal  and  Hanchak,  1999].  The  mate¬ 
rial  parameter  is  described  in  Tabe  1.  Quantitatively, 
the  numerical  result  of  the  present  formulation  cap¬ 
tures  the  plug  velocity  with  good  agreement  as  demon¬ 
strated  in  Figure  1.  The  numerical  model  effectively 
captures  this  plugging  failure  qualitatively  as  shown  in 
Figure  2. 

Concluding  Remarks 


Employing  Cholesky  decomposition  for  the  matrix  Qp, 


Q  p  =  CTVC  (45) 

* 

M  =  (46) 

2°.  If  n  =  0  find  Aqn+i  by 

CTVCAqn+1  =  Atqn  +  ^At2q„  (47) 

3°  Calculate  resultant  of  stresses  within  structure 

rn+i  and  the  vector  rn+i 

Au„+i  =  L-T(Aqn+i)  (48) 

r„+i  =  r(<rra  +  Acr(Aun+i))  (49) 

rn+i  =  L-Trn+i  (50) 

4°  Using  return  map  method  to  get  s'g1 : 

1  =  s«  -  I  ( AgT  +  AAAt)  (51) 


In  this  exposition,  we  described  a  continuum  damage 
mechanics  based  Lagrangian  hydrodynamic  computa¬ 
tional  framework  for  3D  high-  and  hyper-velocity  con¬ 
tact/impact/damage/penetration  applications.  The 
contributions  of  the  paper  included  two  aspects:  (i) 
a  hydrodynamic  constitutive  relation  including  con¬ 
tinuum  damage  mechanics;  (ii)  a  robust  time  inte¬ 
gration  algorithm  for  the  coupled  equation  of  state 
and  the  internal  energy  equation.  A  numerical  ex¬ 
ample  were  performed  to  validate  the  robustness  and 
the  accuracy  of  the  proposed  framework  for  the  con¬ 
tact  / impact / damage / penetration  applications. 
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HY-100  Steel 

4340  Rc  38  steel 

Material  Properties 

Shear  Modulus  [GPa] 

77.5 

77.5 

Density  [kg/m3] 

7801 

7822 

JC  Strength 

A  [MPa] 

758.3 

600 

B  [MPa] 

402.3 

509.5 

C 

0.011 

0.014 

To  [K] 

290 

290 

Tm  [K] 

1817 

1793 

n 

0.26 

0.26 

m 

1.13 

1.03 

cp  [J/kg  K] 

477 

477 

JC  Fracture 

Di 

-0.61 

-0.8 

d2 

2.07 

2.1 

d3 

-0.5 

-0.5 

d4 

0.01 

0.002 

D5 

0 

0.61 

Spall  strength  [GPa] 

5.723 

5.723 

Fracture  strain 

0.035 

0.035 

Damage 

( Pc 

0.99 

0.99 

£d/£F 

0.6 

0.6 

EOS 

Ki  [GPa] 

164 

164 

K2  [GPa] 

294 

294 

K3  [GPa] 

500 

500 

r 

1.16 

1.16 

Table  1:  Material  constants  for  HY-100  Steel  and  4340  Rc  38  steel 


Figure  1:  Velocity  of  plug  node  116445  for  an  impact  velocity  of  246  m/s.  The  nodal  velocity  oscillates  close  to  the 
experimental  value  of  260  m/s. 
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(a)  Forrest al-Hanchak  experimental  results  at  10  /is  intervals 
starting  at  25  /is  for  an  impact  velocity  of  246  m/s  [Forrestal 
and  Hanchak,  1999] 


(b)  t=25  /is 


(c)  t=35  /is 


(d)  t=45  /is  (e)  t=55  /is 


(f)  t=65  /is  (g)  t=75  /is 


(h)  t=85  /is  (i)  t=95  /is 


Figure  2:  Comparison  of  experimental  and  numerical  simulation  results  showing  the  plugging  failure  and  dynamic 
deformation  of  the  HY-100  steel  plate  for  an  impact  velocity  of  246  m/s 
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